#include "TStopwatch.h"

/*
grep MN2D MN2D.h | grep ";" -v | grep class -v  | grep -v grep | awk -F\: '{print $3}'  | sort |wc
grep static MN2D.h | grep return -v | grep "(" | awk '{print $3}'|sort |wc

31
32        //Polint definition
*/

class MN2D
{
  //==============================================================================
  //                                    controls
  //==============================================================================
 public:
  static Bool_t IsMeanFixTrue(){return fkMeanFixTrue;}

 private:
  static Double_t MaxIter(){return 1e8;}

  static Int_t fkParaResolution;
  static Int_t fNPolint; //must be an even number
  static Bool_t fkMeanFixTrue;
  static Bool_t fkResolutionFixTrue;
  static Int_t fRegOpt; 
  //==============================================================================
  //                                    initialization and destruction
  //==============================================================================
 public:
  static void Ini(const Int_t kparares, const Int_t npolint, TH2D * hrawdata, TH1D * hrawyield,  TGraphErrors *greftot, const Int_t nb, const Int_t binx[], const Double_t bbpar[], const Double_t binpar[], const Double_t * respar=0x0);
  static void Clear();
  
 private:
   
  //==============================================================================
  //                                    functions
  //==============================================================================
 public:

 private:
  static Double_t IDdEdxM2Q0(const Double_t *par, const Double_t mean);
  static Double_t IDdEdxM3Q0(const Double_t *par, const Double_t mean);
 
  //==============================================================================
  //                                    parameters
  //==============================================================================
 public:
  static Int_t NparBin(){return dEdxGausN::Npar();}

 private:
  static Int_t NparBB(){return dEdxGausN::NparBB();}//included in Buffer
  static Int_t Npar(){return NparBuffer()+fNbin*NparBin();}
  static Int_t Nfree(){return Npar()-fNfix;}
  static Int_t NparBuffer();
  
  static void IniResPar(Double_t * respar);
  static Int_t NparRes();
  static Double_t GetResolutionParametrization(const Double_t mean, const Double_t * par);
  
  static void SetIfixed();
  static void SetLowHigh();

  static const Double_t * GetBinPar(const Int_t ibin, const Double_t xx, Double_t *localpar, Int_t *outorder, const Bool_t kprint=0);
  static void SetPhysics(const TMatrixDSym * fullcov);

  static Int_t fNfix;
  static Bool_t *fIfixed;
  static Double_t *fLows;
  static Double_t *fHighs;
  static Double_t *fPars;
  static Double_t *fLowErrs;
  static Double_t *fHighErrs;
  static Double_t *fAs;
  static Double_t *fErrAs;
  static Double_t *fRatio;
  static Double_t *fErrRatio;
  static Double_t *fACov;
  //==============================================================================
  //                                    output
  //==============================================================================
 public:
  static void Draw2d(const TString tag);

 private:
  static Char_t* Status(Char_t * info, const Bool_t kend=kTRUE);
  static void Print(const TString tag, const Bool_t kprintBinPar=kTRUE, const Bool_t kprintBinErr=kTRUE, const Bool_t ktofile=kTRUE);
  static void PrintPhysics(FILE *pout, const TString ss, const Double_t * var);
  static void Save(const TString tag);
  static void SavePhysics(TTreeSRedirector * stream, const TString tn, const Int_t ibin);
  static FILE * fLog;

  //==============================================================================
  //                                    fit
  //==============================================================================
 public:
  //(old) Fit and VirtualFit removed from r218 on.
  //MathMinimizer removerd from r308 on.
  //OLDfit removed from r681 on
  static void Fit(const Int_t startnbeta=0);

 private:
  static Double_t LinearFit(const Int_t np, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget, const Int_t order);
  //static Double_t GetSmoothEvalY(const Int_t npo, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget);
  static Double_t GetXX(const Int_t ix);
  static Double_t GetFCNChisquare(Double_t pars[]);
  static void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  static Int_t DoMigrad(const Int_t nbeta, Int_t & counter, const Double_t tol, const Int_t stra);
  static void BetaFit(Int_t & nbeta);
 
  static Int_t NDF(){return fNdata - Nfree();}
  static Double_t ReducedChi2(){ if(!NDF()) exit(1); return fStatChi2/NDF();}
  static Double_t myFCN(){ return fStatChi2 + fBeta*fPenalty; }

  //static TLinearFitter *fFitter;
  static TH2D *fHFitData;
  static TH1D *fHFitYield;
  static TH1D *fHRawYield;
  static TGraphErrors *fGRefTot;

  static Int_t fNbin;
  static Int_t *fBinX;
  static Int_t fNdata;

  static Double_t fStatChi2;
  static Double_t fBeta;
  static Double_t fPenalty;

  static TStopwatch * fStopwatch;
  static TMinuit *fMinuit;
  static Int_t fkErr;

  //==============================================================================
  //                                    end
  //==============================================================================
};

Int_t MN2D::fkParaResolution = 0;
Int_t MN2D::fNPolint = -999;

Bool_t MN2D::fkMeanFixTrue = kFALSE;
Bool_t MN2D::fkResolutionFixTrue = kFALSE;
Int_t MN2D::fRegOpt = 0;

//TLinearFitter * MN2D::fFitter = 0x0;
TH2D *MN2D::fHFitData = 0x0;
TH1D *MN2D::fHFitYield=0x0;
TH1D *MN2D::fHRawYield=0x0;
TGraphErrors *MN2D::fGRefTot=0x0;

Int_t MN2D::fNbin = 0;
Int_t *MN2D::fBinX = 0x0;
Int_t MN2D::fNdata = 0;

Int_t MN2D::fNfix = 0;
Bool_t *MN2D::fIfixed=0x0;
Double_t *MN2D::fLows=0x0;
Double_t *MN2D::fHighs=0x0;
Double_t *MN2D::fPars=0x0;
Double_t *MN2D::fLowErrs=0x0;
Double_t *MN2D::fHighErrs=0x0;
Double_t *MN2D::fAs=0x0;
Double_t *MN2D::fErrAs=0x0;
Double_t *MN2D::fRatio=0x0;
Double_t *MN2D::fErrRatio=0x0;
Double_t *MN2D::fACov=0x0;

Double_t MN2D::fStatChi2=-999;
Double_t MN2D::fBeta = -999;
Double_t MN2D::fPenalty=-999;

FILE * MN2D::fLog = 0x0;
TStopwatch *MN2D::fStopwatch = 0x0;
TMinuit *MN2D::fMinuit = 0x0;
Int_t MN2D::fkErr = -999;

//==============================================================================
//                                    initialization and destruction
//==============================================================================

void MN2D::Ini(const Int_t kparares, const Int_t npolint, TH2D * hrawdata, TH1D * hrawyield, TGraphErrors *greftot, const Int_t nb, const Int_t binx[], const Double_t bbpar[], const Double_t binpar[], const Double_t * respar)
{
  fkParaResolution = kparares;

  const TString pwd = gSystem->pwd();

  if(npolint%2){
    printf("MN2D::Ini npolint must be even!! %d", npolint);exit(1);
  }
  fNPolint = npolint;

  fHRawYield = (TH1D*) hrawyield->Clone("hrawyield");
  dEdxGausN::GetHFit(hrawdata, hrawyield, fHFitData, fHFitYield);

  if(greftot){
    fGRefTot = (TGraphErrors*) greftot->Clone("mn2dreftot");
  }
  if(dEdxGausN::IsWithRef() && !fGRefTot){
    printf("MN2D::Int fGRefTot null!! %d\n", dEdxGausN::IsWithRef()); exit(1);
  }

  fNbin = nb;

  fBinX = new Int_t[fNbin];
  memcpy(fBinX, binx, sizeof(Int_t)*fNbin); 
  fHFitData->GetXaxis()->SetRange(fBinX[0], fBinX[fNbin-1]);

  fIfixed = new Bool_t[Npar()];

  if(pwd.Contains("RegOpt0")){
    fRegOpt=0;
  }
  else if(pwd.Contains("RegOpt1")){
    fRegOpt=1;
  }
  else if(pwd.Contains("RegOpt2")){
    fRegOpt=2;
  }
  else if(pwd.Contains("RegOpt3")){
    fRegOpt=3;
  }
  else if(pwd.Contains("RegOpt4")){
    fRegOpt=4;
  }
  else{
    printf("MN2D::Ini default RegOpt %d\n", fRegOpt);
  } 

  if(pwd.Contains("MeanFixTrue")){
    fkMeanFixTrue = kTRUE;
  }
  else{
    fkMeanFixTrue = kFALSE;
  }

  if(dEdxGausN::IsWithRef()){
    fkMeanFixTrue = kTRUE;
  }

  if(pwd.Contains("ResolutionFixTrue")){
    fkResolutionFixTrue = kTRUE;
  }
  else{
    fkResolutionFixTrue = kFALSE;
  }
  SetIfixed();

  fLows = new Double_t[Npar()];
  fHighs= new Double_t[Npar()];
  SetLowHigh();

  fPars = new Double_t[Npar()];
  memcpy(fPars,              bbpar,  sizeof(Double_t)*NparBB());

  if(respar){
    Double_t * tmppar = &(fPars[NparBB()]);
    tmppar[0]=respar[0];
    tmppar[1]=respar[1];
    tmppar[2]=respar[2];
  }
  else{
    IniResPar(&(fPars[NparBB()]));
  }

  memcpy(&(fPars[NparBuffer()]), binpar, sizeof(Double_t)*fNbin*NparBin());

  fLowErrs = new Double_t[Npar()];
  fHighErrs = new Double_t[Npar()];

  for(Int_t ii=0; ii<Npar(); ii++){
    fLowErrs[ii]=0;
    fHighErrs[ii]=0;
  }

  //only correlation of A's with pion is saved!!
  fAs = new Double_t[fNbin*dEdxGausN::Ntype()];
  fErrAs = new Double_t[fNbin*dEdxGausN::Ntype()];
  fRatio = new Double_t[fNbin*dEdxGausN::Ntype()];
  fErrRatio = new Double_t[fNbin*dEdxGausN::Ntype()];
  fACov = new Double_t[fNbin*dEdxGausN::Ntype()];

  for(Int_t ii=0; ii<fNbin*dEdxGausN::Ntype(); ii++){
    fAs[ii]=0;
    fErrAs[ii]=0;
    fRatio[ii]=0;
    fErrRatio[ii]=0;
    fACov[ii]=0;
  }

  fLog = fopen("logmn2d.txt","w");
  fStopwatch = new TStopwatch;

  Print("MN2D::Ini stdout pre GetFCNChisquare",1,0,0);

  fBeta = 0.0;
  GetFCNChisquare(fPars);

  Print("MN2D::Ini fLog");
}

void MN2D::Clear()
{
  delete fBinX;
  delete fPars;
  delete fAs;
  delete fErrAs;
  delete fRatio;
  delete fErrRatio;
  delete fACov;
  delete fLowErrs;
  delete fHighErrs;
  delete fIfixed;
  delete fLows;
  delete fHighs;

  fclose(fLog);
  //FILE can not be deleted, otherwise crashes!
  //delete fLog;

  delete fStopwatch;
  delete fMinuit;
}

//==============================================================================
//                                    functions
//==============================================================================
Double_t MN2D::IDdEdxM2Q0(const Double_t *par, const Double_t mean)
{
  //return (TMath::Abs(par[0]) * TMath::Power(mean,par[1]) + par[2] )/mean;

  //rewritting the power term in denominator and add eps to protect!
  return TMath::Abs(par[0])/( TMath::Power(mean,1-par[1])+EPSILON )+ par[2]/mean;
}

Double_t MN2D::IDdEdxM3Q0(const Double_t *par, const Double_t mean)
{
  return TMath::Abs(par[0])/TMath::Sqrt(mean) + par[1] + par[2]/mean;
}

//==============================================================================
//                                    parameters
//==============================================================================
Int_t MN2D::NparBuffer()
{
  return NparBB()+NparRes();
}

Int_t MN2D::NparRes()
{
  if(fkParaResolution%100==0){//{2,3}{0,1}00
    return 3;
  }
  else{
    printf("MN2D::NparRes abnormal fkParaResolution!! %d\n", fkParaResolution);exit(1);
  }
}


void MN2D::SetIfixed()
{
  for(Int_t ipar=0; ipar<Npar(); ipar++){
    fIfixed[ipar] = kFALSE;
  }
  fNfix = 0;

  if(fkMeanFixTrue){
    for(Int_t ii=0; ii<NparBB(); ii++){
      fIfixed[ii]=kTRUE; fNfix++;
    }
  }

  //for fkParaResolution = 0,  NparRes()=0
  if(fkResolutionFixTrue){
    for(Int_t ip=0; ip<NparRes(); ip++){
      fIfixed[NparBB()+ip] = kTRUE; fNfix++;
    }
  }

  for(Int_t ibin=0; ibin<fNbin; ibin++){
    Bool_t *mgpfix=&(fIfixed[NparBuffer()+ibin*NparBin()]);

    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      mgpfix[itype*dEdxGausN::NparType()+1] = kTRUE; fNfix++;
      mgpfix[itype*dEdxGausN::NparType()+2] = kTRUE; fNfix++;
    }
  }
}

void MN2D::SetLowHigh()
{
  for(Int_t ipar=0; ipar<Npar(); ipar++){
    fLows[ipar] = 0;
    fHighs[ipar]= 0;
  }

  Double_t *bslow =&(fLows[ NparBuffer()]);
  Double_t *bshigh=&(fHighs[NparBuffer()]);
  for(Int_t ibin=0; ibin<fNbin; ibin++){
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      //For variable parameters with limits, MINUIT uses the following transformation:
      //P   = arcsin(2((P   -a)/(b- a))-1)                P    = a+((b- a)/(2))(sinP    +1)

      //limits imposed in fith already, so DON'T use them here
      //---> important to have limits, otherwise kaon stay at 0
      //---> limits will be released in nbeta>0
      bslow[ ibin*NparBin()+itype*dEdxGausN::NparType()+0] = 0;
      bshigh[ibin*NparBin()+itype*dEdxGausN::NparType()+0] = 1;     
    }
  }
}

void MN2D::IniResPar(Double_t * respar)
{
  Double_t * tmppar = &(respar[0]);

  const Int_t dedxbit = fkParaResolution/1000;
  //dEdx{2,3}xxx
  if(dedxbit==2){
    tmppar[0]=dEdxGausN::IniResolution();
    tmppar[1]=1;
    tmppar[2]=0;
  }
  else if(dedxbit==3){
    tmppar[0]=0;
    tmppar[1]=dEdxGausN::IniResolution();
    tmppar[2]=0;
  }
  else if(dedxbit){
    printf("MN2D::IniResPar wrong fkParaResolution %d\n", fkParaResolution); exit(1);
  }
}

Double_t MN2D::GetResolutionParametrization(const Double_t mean, const Double_t * rawParRes)
{
  //remember to updata NparRes() and IniResPar()

  const Double_t *par = &(rawParRes[0]);

  Double_t res = -999;

  switch(fkParaResolution){
  case 2000:
    res = IDdEdxM2Q0(par, mean);
    break;
  case 3000:
    res = IDdEdxM3Q0(par, mean);
    break;
  default:
    printf("MN2D::GetResolutionParametrization wrong fkParaResolution %d\n", fkParaResolution); exit(1);
  }

  //cut out negative branched to improve minimum searching, even though negative values protected by abs in resolutionGaus
  return TMath::Max(EPSILON, res);
}

const Double_t * MN2D::GetBinPar(const Int_t ibin, const Double_t xx, Double_t *glopar, Int_t *outorder, const Bool_t kprint)
{  
  //
  //glopar is modified 
  //1) in buffer for mean from output after SetMean no matter whenther fixmean=1
  //2) and resolution from GetResolutionParameerization

  if(kprint){
    printf("getbinpar before %e %e %e %e %e\n", glopar[0], glopar[1], glopar[2], glopar[3], glopar[4]);
  }
  Double_t *binpar=&(glopar[NparBuffer()+ibin*NparBin()]);

  Int_t myorder[dEdxGausN::Ntype()];
 
  //============= from buffer BB part to ibin Mean in binpar
  //put bbmean to binpar
  //put order to myorder
  //glopar (BBpar) not changed

  if(kprint){
    printf("getbinpar after %e %e %e %e %e\n", glopar[0], glopar[1], glopar[2], glopar[3], glopar[4]);
  }
  dEdxGausN::SetMean(xx, glopar, binpar, myorder, kprint);

  //============= resolution from buffer just after BB
  Double_t * respar = &(glopar[NparBB()]);

  //=============
  Double_t realmean = 0;
  for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
    const Int_t gausanchor = ith*dEdxGausN::NparType();

    //>>>>>>>>>>>>>>>>> calc. width and set width to glopar(i.e. binpar) <<<<<<<<<<<<<<<<<<<<<<<
    //big bug fixed!!
    realmean += binpar[gausanchor+1];
    
    //input is itype, i.e. myorder[ith]=itype
    const Double_t calcRes = GetResolutionParametrization(realmean, respar);
    binpar[gausanchor+2] = calcRes;    
  }

  if(outorder){
    memcpy(outorder,myorder, sizeof(myorder));
  }

  if(kprint){
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      for(Int_t ip=0; ip<dEdxGausN::NparType(); ip++){
        const Int_t ipar = ith*dEdxGausN::NparType()+ip;
        printf("ibin %d xx %e ith %d ip %d ipar %d -- %e\n", ibin, xx, ith, ip, ipar, binpar[ipar]);
      }
    }
  }
  return binpar;
}

void MN2D::SetPhysics(const TMatrixDSym * fullcov)
{
  //in each pbin, do
  //1. get the order, find where is pion
  //2. set As and ErrAs, should be identical to the one in ::Save
  //3. calc. ratio over pion
  //4. calc. err of ratio from fullcov
  //5. calc. sumcov

  for(Int_t ibin = 0; ibin < fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = GetXX(ix);

    Double_t dummy[NparBin()];//big bug fixed!!
    //itype = order[ith]
    Int_t order[dEdxGausN::Ntype()];
    dEdxGausN::SetMean(xx, fPars, dummy, order);

    //ith=invmap[itype]
    Int_t invmap[dEdxGausN::Ntype()];
    dEdxGausN::GetInvMap(order, invmap);

    Double_t sumcov=0;
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      const Int_t isave = ibin*dEdxGausN::Ntype()+order[ith];
      const Int_t ipar  = NparBuffer() + ibin*NparBin() + ith*dEdxGausN::NparType();
      //printf("longtest %d %d -- %d %d -- %d\n", ith, order[ith], ipar, ibin*dEdxGausN::Ntype()+order[ith], ibin);

      fAs[isave] = TMath::Abs(fPars[ipar]);
     
      const Double_t sigma2 = (*fullcov)[ipar][ipar];
      if(sigma2<EPSILON){
        if(TMath::Abs(sigma2)<EPSILON){
          fErrAs[isave] = 0;
        }
        else{
          printf("MN2D::SetPhysics sigma2<EPSILON %e %d\n", sigma2, ipar); exit(1);
        }
      }
      else{
        fErrAs[isave] = TMath::Sqrt(sigma2);
      }

      const Int_t ipion = NparBuffer() + ibin*NparBin() + invmap[dEdxGausN::kPION]*dEdxGausN::NparType();
      fACov[isave] = (*fullcov)[ipar][ipion];

      //sum of cov elements and then sqrt is error fo \sum A
      for(Int_t jth=0; jth<dEdxGausN::Ntype(); jth++){
        const Int_t jpar   = NparBuffer() + ibin*NparBin() + jth*dEdxGausN::NparType();
        sumcov+=(*fullcov)[ipar][jpar];
      }
    }

    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      const Int_t ia = ibin*dEdxGausN::Ntype()+itype;
      const Int_t ib = ibin*dEdxGausN::Ntype()+dEdxGausN::kPION;
      const Double_t aa = fAs[ia];
      const Double_t bb = fAs[ib];
      if(bb<EPSILON){
        fRatio[ia]=fErrRatio[ia]=-999;
      }
      else{
        const Double_t ff = aa/bb;
        const Double_t v0 = ff/aa;
        const Double_t v1 = -ff/bb;
        fRatio[ia]=ff;
        fErrRatio[ia]=BaseUtils::GetCorrelatedError(v0,v1,TMath::Power(fErrAs[ia],2), fACov[ia], TMath::Power(fErrAs[ib],2));
      }
    }
  }
}

//==============================================================================
//                                    output
//==============================================================================
Char_t* MN2D::Status(Char_t * info, const Bool_t kend)
{
  Char_t * status = Form("fkErr = %d, fStatChi2 = %e, NDF() = %d, ReducedChi2() = %e, fBeta =  %e, fPenalty = %e, b*p = %e fcn = %e", fkErr, fStatChi2, NDF(), ReducedChi2(), fBeta, fPenalty, fBeta*fPenalty, myFCN());
  printf("%-110s %s%c", info, status, kend?'\n':' ');
  return status;
}

void MN2D::Print(const TString tag, const Bool_t kprintBinPar, const Bool_t kprintBinErr, const Bool_t ktofile)
{
  FILE *pout = stdout;
  if(ktofile)
    pout=fLog;

  fprintf(pout, "\n================================================ %s:\n", tag.Data());

  fprintf(pout, "fkParaResolution = %d;\n", fkParaResolution);
  fprintf(pout, "fNPolint = %d;\n", fNPolint);
  fprintf(pout, "fRegOpt = %d;\n", fRegOpt);
  fprintf(pout, "fkMeanFixTrue = %d;\n", fkMeanFixTrue);
  fprintf(pout, "fkResolutionFixTrue = %d;\n", fkResolutionFixTrue);
  
  fprintf(pout, "NparBuffer(): %d;\n", NparBuffer());
  fprintf(pout, "NparBB(): %d;\n", NparBB());
  fprintf(pout, "NparRes(): %d\n", NparRes());

  fprintf(pout, "NparBin(): %d;\n", NparBin());
  fprintf(pout, "Npar(): %d;\n", Npar());
  fprintf(pout, "Nfree(): %d\n", Nfree());
  fprintf(pout, "fNfix = %d\n", fNfix);
  fprintf(pout, "fNdata = %d\n", fNdata);
  fprintf(pout, "Status: %s\n", Status(Form("-> MN2D::Print %s", tag.Data())));

  fprintf(pout, "fNbin = %d;\n", fNbin);

  for(Int_t ii=0; ii<NparBuffer(); ii++){
    fprintf(pout, "fPars[%d] = %e;\n", ii, fPars[ii]);
  }
  
  fprintf(pout, "\n");

  if(!kprintBinPar)
    return;

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fPars[%3d] = %10e; ", id, fPars[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  if(!kprintBinErr)
    return;

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fLowErrs[%3d] = %10e; ", id, fLowErrs[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  for(Int_t ii=0; ii<fNbin; ii++){
    fprintf(pout, "fBinX[%3d] = %d; ", ii, fBinX[ii]);
    for(Int_t ip = 0; ip<NparBin(); ip++){
      const Int_t id = NparBuffer()+ii*NparBin()+ip;
      fprintf(pout, "fHighErrs[%3d] = %10e; ", id, fHighErrs[id]);
    }
    fprintf(pout, "\n");
  }

  fprintf(pout, "\n");

  PrintPhysics(pout, "fAs", fAs);
  PrintPhysics(pout, "fErrAs", fErrAs);
  PrintPhysics(pout, "fACov", fACov);
  PrintPhysics(pout, "fRatio", fRatio);
  PrintPhysics(pout, "fErrRatio", fErrRatio);  

  fprintf(pout, "\n");
}

void MN2D::PrintPhysics(FILE *pout, const TString ss, const Double_t * var)
{
  for(Int_t ibin=0; ibin<fNbin; ibin++){
    fprintf(pout, "bin=%4d; ", ibin);
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      const Int_t id = ibin*dEdxGausN::Ntype()+itype;
      fprintf(pout, "%s[%3d] = %10e; ", ss.Data(), id, var[id]);
    }
    fprintf(pout, "\n");
  }
  fprintf(pout, "\n");
}

void MN2D::Save(const TString tag)
{
  TTreeSRedirector *fStream = new TTreeSRedirector(Form("%s_save_%s.root", tag.Data(), dEdxGausN::Opt()));

  for(Int_t ip=0; ip<NparBuffer(); ip++){
    (*fStream)<<"tree"<<
      Form("bb%d=",ip)<<fPars[ip]<<
      Form("leb%d=",ip)<<fLowErrs[ip]<<
      Form("heb%d=",ip)<<fHighErrs[ip];
  }
  Double_t chidof = ReducedChi2();
  Int_t npar = Npar();
  Int_t ndf = NDF();
  (*fStream)<<"tree"<<
    "npar="<<npar<<
    "chi2="<<fStatChi2<<
    "ndof="<<ndf<<
    "chidof="<<chidof<<
    "penalty="<<fPenalty<<
    "nbin="<<fNbin<<
    "\n";

  //localpar, i.e. fPars will be updated in GetBinPar (bufferResPart, mean and resolution are set).
  Double_t *localpar = fPars;
  Double_t *locallowerr  = fLowErrs;
  Double_t *localhigherr = fHighErrs;

  for(Int_t ibin=0; ibin<fNbin; ibin++){
    Int_t ix = fBinX[ibin];
    Double_t xx = GetXX(ix);

    //need to put order here, otherwise the output has wrong order
    //only record order, no reordering
    Int_t order[dEdxGausN::Ntype()];
    GetBinPar(ibin, xx, localpar, order);

    //totalyield, outpar, outlowerr, outhigherr have to be defined here due to streamer
    Double_t yields[dEdxGausN::Ntype()];
    Double_t outpar[npar];
    Double_t outlowerr[npar];
    Double_t outhigherr[npar];
    Double_t totalyield = fHRawYield->GetBinContent(ix);
    const Int_t ipar = NparBuffer()+ibin*NparBin();
    dEdxGausN::SavePar(fStream, "mn2d", &totalyield, yields, &(localpar[ipar]), &(locallowerr[ipar]), &(localhigherr[ipar]), order, outpar, outlowerr, outhigherr);

    SavePhysics(fStream,"mn2d", ibin);

    (*fStream)<<"mn2d"<<
      "ix="<<ix<<
      "xx="<<xx<<   
      "\n";
  }
  
  fStream->GetFile()->cd();
  fHFitData->Write();
  delete fStream;

  Print(Form("MN2D::Save %s", tag.Data()));
}

void MN2D::SavePhysics(TTreeSRedirector * stream, const TString tn, const Int_t ibin)
{
  for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
    const Int_t id = ibin*dEdxGausN::Ntype()+itype;
    (*stream)<<tn<<
      Form("aa%d=",   itype)<<fAs[id]<<
      Form("ea%d=",   itype)<<fErrAs[id]<<
      Form("cova%d=", itype)<<fACov[id]<<
      Form("rr%d=",   itype)<<fRatio[id]<<
      Form("er%d=",   itype)<<fErrRatio[id];
  }
}

void MN2D::Draw2d(const TString tag)
{
  //_______________________________________ define output _____________________
  TList *ll = new TList;
  ll->SetOwner(kTRUE);

  TH2D *horiginal = (TH2D*)fHFitData->Clone(tag+"_original");
  ll->Add(horiginal);

  TH2D *hfit = (TH2D*)fHFitData->Clone(tag+"_fit");
  hfit->Scale(0);
  ll->Add(hfit);

  const Double_t * fitxarr= BaseUtils::GetAxisArray(hfit->GetXaxis());
  if(!fitxarr){
    printf("MN2D::Draw2d null fitxarr %d!!\n", hfit->GetNbinsX()); exit(1);
  }

  TH1D *hmean[dEdxGausN::Ntype()];
  for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
    hmean[itype] = new TH1D(tag+Form("_mean%d", itype),"",hfit->GetNbinsX(), hfit->GetXaxis()->GetXmin(), hfit->GetXaxis()->GetXmax());
    //big bug fixed!
    hmean[itype]->GetXaxis()->Set(hfit->GetNbinsX(), fitxarr);
 
    //hmean[itype]->GetXaxis()->SetRange(fBinX[0], fBinX[fNbin-1]);
    ll->Add(hmean[itype]);
  }

  TH1D * hchi2 = new TH1D(tag+"_chi2","",hfit->GetNbinsX(), hfit->GetXaxis()->GetXmin(), hfit->GetXaxis()->GetXmax());
  //big bug fixed!
  hchi2->GetXaxis()->Set(hfit->GetNbinsX(), fitxarr);

  hchi2->GetXaxis()->SetRange(fBinX[0], fBinX[fNbin-1]);
  hchi2->SetMarkerStyle(20);
  ll->Add(hchi2);

  //______________________________________________________________________
  //______________________________________________________________________

  //localpar is updated in GetBinPar, (bufferResPart, mean and resolution is set).
  Double_t *localpar = fPars;

  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = GetXX(ix);

    Int_t order[100];
    const Double_t * binpar = GetBinPar(ibin, xx, localpar, order);

    //resotre to natural gaus par, order not changed!!
    Double_t gauspar[NparBin()];
    dEdxGausN::ToGausPar(binpar, gauspar);
    Int_t invmap[100];
    dEdxGausN::GetInvMap(order, invmap);
    
    Double_t chi2bini = 0;
    Double_t ndfbini = 0;
    for(Int_t iy = 1; iy<= hfit->GetNbinsY(); iy++){
      const Double_t yy = hfit->GetYaxis()->GetBinCenter(iy);
      const Double_t cc = dEdxGausN::Func(&yy, binpar);

      //set hfit
      hfit->SetBinContent(ix, iy, cc);

      //calc chi2
      const Double_t rawc = fHFitData->GetBinContent(ix, iy);
      const Double_t rawerr = fHFitData->GetBinError(ix, iy);
      if(rawerr>EPSILON){
        chi2bini+=TMath::Power((rawc-cc)/rawerr,2);
        ndfbini++;
      }

      //set mean
      for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
        hmean[itype]->SetBinContent(ix, gauspar[invmap[itype]*dEdxGausN::NparType()+1]);
        hmean[itype]->SetBinError(ix, 0);
      }
    }

    //set chi2
    ndfbini -= Nfree()/fNbin  ;
    hchi2->SetBinContent(ix, ndfbini>EPSILON? chi2bini/ndfbini : 0);
    hchi2->SetBinError(ix, 0);

    //get slice
    TH1D *tmph = fHFitData->ProjectionY(Form("%sx%d", fHFitData->GetName(), ix),ix,ix,"oe");
    tmph->SetTitle(Form("x %e",xx));
    tmph->SetLineColor(kBlack);
    tmph->SetMarkerColor(kBlack);
    ll->Add(tmph);

    //draw slice
    TH1D *hslicefit = BaseUtils::GetHfit(Form("%sfit", tmph->GetName()), dEdxGausN::Func, binpar, tmph->GetXaxis()->GetXmin(), tmph->GetXaxis()->GetXmax());
    hslicefit->SetLineColor(kBlack);
    //hslicefit->SetDrawOption("same hist C");
    ll->Add(hslicefit);

    Int_t col[]={kBlue, kRed, kGreen+3, kMagenta, kYellow};
    for(Int_t itype=0; itype<dEdxGausN::Ntype(); itype++){
      dEdxGausN::DrawTF1(itype, tmph->GetName(), &(gauspar[invmap[itype]*dEdxGausN::NparType()]), col[itype], ll, fLog); 
    }
  }

  //remove empty bins
  TGraph * grmean[dEdxGausN::Ntype()];
  for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
    grmean[itype] = new TGraph;
    grmean[itype]->SetName(Form("%sgr",hmean[itype]->GetName()));
    ll->Add(grmean[itype]);
    Int_t ig = 0;
    for(Int_t ii=1; ii<=hmean[itype]->GetNbinsX(); ii++){
      const Double_t yy = hmean[itype]->GetBinContent(ii);
      if(yy<EPSILON){
        continue;
      }

      const Double_t xx = hmean[itype]->GetBinCenter(ii);

      grmean[itype]->SetPoint(ig, xx, yy);
      ig++;
    }
  }

  TGraph *grchi2 = new TGraph;
  grchi2->SetName(Form("%sgr",hchi2->GetName()));
  ll->Add(grchi2);
  Int_t ig = 0;
  for(Int_t ii=1; ii<=hchi2->GetNbinsX(); ii++){
    const Double_t yy = hchi2->GetBinContent(ii);
    if(yy<EPSILON){
      continue;
    }
    
    const Double_t xx = hchi2->GetBinCenter(ii);
    
    grchi2->SetPoint(ig, xx, yy);
    ig++;
  }

  TFile *fout= new TFile(Form("%s_draw2d_%s.root", tag.Data(), dEdxGausN::Opt()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
  delete fout;
  delete ll;

  Print(Form("MN2D::Draw2d %s", tag.Data()));
}

//==============================================================================
//                                    fit
//==============================================================================
/*
Double_t MN2D::GetSmoothEvalY(const Int_t npo, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget)
{
  if(!fFitter){
    fFitter = new TLinearFitter(1,"1++x++x*x");
  }

  fFitter->ClearPoints();
  for(Int_t ii=0; ii<npo; ii++){
    Double_t xi = xarr[ii]-xtarget;
    fFitter->AddPoint(&xi, yarr[ii], eys[ii]<EPSILON ? EPSILON : eys[ii]);
  }
  fFitter->EvalRobust();
  return fFitter->GetParameter(0);
}
*/

Double_t MN2D::LinearFit(const Int_t np, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget, const Int_t order)
{
  //chi2 polinomial fit of degree, i.e. up to x^(deg-1), see also Recipes p673.
  const Int_t deg = order+1;
  TMatrixD leftmatrix(deg, deg);
  TMatrixD rightvector(deg,1);

  //loop over data point
  for(Int_t ip=0; ip<np; ip++){
    const Double_t xx = xarr[ip]-xtarget;
    const Double_t yy = yarr[ip];
    //to avoid ey<0 or ey=0
    const Double_t ey = TMath::Max(EPSILON, eys[ip]);

    //starting from [0][0] term
    Double_t lastrecord = 1/ey/ey;    
    //loop over each matrix entry once with increasing power of xx
    for(Int_t sumij=0; sumij<=(deg-1)*2; sumij++){
      for(Int_t ii=TMath::Max( 0, sumij-(deg-1) ); ii<=TMath::Min(deg-1, sumij); ii++){
        const Int_t jj = sumij - ii;
        //printf("%d: %d %d\n", sumij, ii, jj);
        
        leftmatrix[ii][jj] += lastrecord;
      }
      //printf("sum %d\n\n",sumij);
      lastrecord *= xx;
    }

    //loop over vector entry
    lastrecord = yy/ey/ey;
    for(Int_t ii=0; ii<deg; ii++){
      rightvector[ii][0] +=  lastrecord;
      lastrecord *= xx;
    }

    //checked by hand, correct!
    //leftmatrix.Print();
    //printf("===\n");
    //rightvector.Print();
  }

  /*
TMatrixD a(2,2)
a[0][0]=cos(3)
a[0][1]=sin(3)
a[1][0]=-sin(3)
a[1][1]=cos(3)
a.Print()
TDecompLU lu(a)
TMatrixD c=lu.Invert()
c.Print()
a.Print()
  */

  //leftmatrix
  TDecompSVD mdec(leftmatrix);
  if(!mdec.Decompose()){
    printf("MN2D::LinearFit matrix singular?\n");
    for(Int_t ii=0; ii<np; ii++){
      printf("%d: %e %e %e - %e\n", ii, xarr[ii], yarr[ii], eys[ii], xtarget);
    }
    leftmatrix.Print();
    exit(1);
  }
  const TMatrixD inv=mdec.Invert();

  Double_t par0 = 0;
  for(Int_t ii=0; ii<deg; ii++){
    par0 += inv[0][ii]*rightvector[ii][0];
  }

  //compared to LinearFitter::Eval, consistent
  return par0;
}

Double_t MN2D::GetXX(const Int_t ix)
{
  return fHFitData->GetXaxis()->GetBinCenter(ix);
}

Double_t MN2D::GetFCNChisquare(Double_t localpar[])
{
  static Int_t istep = 0;
  const Int_t kprint = 0;//(istep>= 112000 && istep < 113000);// 0;//TMath::Max(istep-46649,0);// 
  if(kprint) {
    for(int ii=0; ii<Npar(); ii++){
      printf("istep %d %d: %e %e\n", istep, ii, fPars[ii], localpar[ii]);
    }
  }

  //statistical chi2, i.e. the standard one
  fStatChi2 = 0;

  fNdata = 0;

  //record fraction for particle types and position
  Double_t frac[dEdxGausN::Ntype()][fNbin];
  Double_t px[fNbin];

  for(Int_t ibin=0;ibin<fNbin; ibin++){
    const Int_t ix = fBinX[ibin];
    const Double_t xx = GetXX(ix);

    Int_t outorder[dEdxGausN::Ntype()];
    //order[ith]=itype
    const Double_t *binpar = GetBinPar(ibin, xx, localpar, outorder, kprint);

    //calc chi2
    const Double_t norm = fHFitYield->GetBinContent(ix);
    for(Int_t iy = 1; iy<= fHFitData->GetNbinsY(); iy++){
      const Double_t yy = fHFitData->GetYaxis()->GetBinCenter(iy);
      //no change of binpar in ::Func
      const Double_t zfit = dEdxGausN::Func(&yy, binpar);

      const Double_t zdata = fHFitData->GetBinContent(ix, iy);

      const Double_t prob = TMath::Poisson(zdata*norm, zfit*norm);
      fStatChi2 += -TMath::Log(prob>EPSILON ? prob : EPSILON);

      fNdata ++;
    }

    if(kprint){
      printf("istep %d binpar %e %e %e %e %e %e outorder %d %d %d %d\n", istep, binpar[0],binpar[1],binpar[2],binpar[3],binpar[4],binpar[5],outorder[0],outorder[1],outorder[2],outorder[3]);
    }

    //save frac
    for(Int_t ith=0; ith<dEdxGausN::Ntype(); ith++){
      //big bug fixed !! only abs() matters
      frac[outorder[ith]][ibin]=TMath::Abs(binpar[ith*dEdxGausN::NparType()+0]);
      //if(kprint) {
      //  printf("%d frac %d %d: %d %e %e\n", istep, ibin, ith, outorder[ith], frac[outorder[ith]][ibin], binpar[ith*dEdxGausN::NparType()+0]);
      //}
    }
    //save px
    px[ibin]=xx;
  }

  //-----------penaltiy term on fractions ----------
  //regularized on all particles, otherwise bad behavior of the free one
 
  fPenalty = 0;

  //fNPolint = 0 no penalty, checked: consistent with before
  //need to be calculated even if beta=0 in order to evaluate next beta
  if(fNPolint){
    //for NPolint = 2
    //* X *
    //  * X *

    //for NPolint = 4
    //* * X * *, 4-point interpolation
    //0, 1, 2, ..., n-3, n-2, n-1
    //    * * X * *

    //ibin must increase with step fNpolint/2, like NPolint = 2
    //otherwise to strong overlap and therefore correlation between neighbouring bins -> never converge
    //----> it seems fine to ues ibin++
    for(Int_t ibin = fNPolint/2; ibin <= fNbin-1-fNPolint/2; ibin++){
      const Double_t targetX = px[ibin];
      Double_t Xpol[fNPolint];
      for(Int_t ip=0; ip<fNPolint/2; ip++){
        ///think about ip=0
        Xpol[ip]           = px[ibin-fNPolint/2+ip];
        Xpol[fNPolint-1-ip]= px[ibin+fNPolint/2-ip];
      }

      const Int_t ix = fBinX[ibin];
      const Double_t reduYield = fHFitYield->GetBinContent(ix);

      for(Int_t itype = 0; itype<dEdxGausN::Ntype(); itype++){
        const Double_t targetY = frac[itype][ibin];
       
        Double_t Ypol[fNPolint];
        Double_t EyPol[fNPolint];
        for(Int_t ip=0; ip<fNPolint/2; ip++){
          Ypol[ip]            = frac[itype][ibin-fNPolint/2+ip];
          Ypol[fNPolint-1-ip] = frac[itype][ibin+fNPolint/2-ip];
          //if(kprint) {
          //  printf("longtest1 %d -- %d %d %d %d %d -- %e %e %e %e\n", itype, fNPolint, ip, ibin-fNPolint/2+ip, fNPolint-1-ip, ibin+fNPolint/2-ip, Ypol[ip]  , Ypol[fNPolint-1-ip], frac[itype][ibin+fNPolint/2-ip],  frac[itype][ibin-fNPolint/2+ip]);
          //}
        }

        for(Int_t ip=0; ip<fNPolint; ip++){
          EyPol[ip] = TMath::Sqrt(Ypol[ip]/reduYield);
        }

        //it should be very fast

        //static void Polint(const Int_t n, const Double_t *xa, const Double_t *ya, const Double_t x, Double_t &y, Double_t &dy, const Int_t kprint=0);

        Double_t Yeva = -999;
        if(fRegOpt==0){
          Double_t dummy = -999;
          BaseUtils::Polint(fNPolint, Xpol, Ypol, targetX, Yeva, dummy);
        }
        else{
          //use pol2 instead, based on Taylor expansion
          Yeva = LinearFit(fNPolint, Xpol, Ypol, EyPol, targetX, fRegOpt);
          if(kprint){
            for(Int_t ipr=0; ipr<fNPolint; ipr++){
              printf("istep %d: %d %e %e %e -- %e %e\n", istep, ipr, Xpol[ipr], Ypol[ipr], EyPol[ipr], targetX, Yeva);
            }
          }
        }

        //chi^2 form of penalty
        const Double_t delta = Yeva-targetY;
        //the Abs() is important, otherwise there is "nan" and not converge
        //use (targetY+Yeva)/2 instead of targetY for <frac>
        const Double_t denom   = TMath::Sqrt(TMath::Abs(targetY+Yeva)/2/reduYield); //= sqrt(AN)/N

        if(denom<EPSILON*1e-6){
          for(Int_t ibad=0; ibad<fNbin; ibad++){
            printf("bad %d %e -- %e %e %e %e\n", ibad, px[ibad], frac[0][ibad], frac[1][ibad], frac[2][ibad], frac[3][ibad]);
          }
          printf("MN2D::GetFCNChisquare denom 0 !! ib %d itype %d tX %e X %e %e %e %e itype %d Y %e %e %e %e -- tY %e Yeva %e delta %e denom %e fPenalty %e\n", ibin, itype, targetX, Xpol[0], Xpol[1], Xpol[2], Xpol[3], itype, Ypol[0], Ypol[1], Ypol[2], Ypol[3], targetY, Yeva, delta, denom, fPenalty); //exit(1);
        }
        else{
          //1/2 to restore to Gaussian p.d.f, note the sigma is considered as parameter-independent
          fPenalty += TMath::Power(delta/denom, 2) / 2. ;
        }
        //poisson form of penalty
        /*
        const Double_t penData = targetY*reduYield;
        const Double_t penFit = Yeva*reduYield;
        const Double_t penP = TMath::Poisson(penData<0?0 : penData, penFit<0?0: penFit );
        fPenalty += -TMath::Log(penP>EPSILON ? penP : EPSILON);
        */

        //if(kprint ){
        //  printf("longtest %d ib %d tX %e X %e %e %e %e itype %d Y %e %e %e %e -- tY %e Yeva %e delta %e denom %e del/den %e w2 %e fPenalty %e -- redy %e\n", istep,  ibin, targetX, Xpol[0], Xpol[1], Xpol[2], Xpol[3], itype, Ypol[0], Ypol[1], Ypol[2], Ypol[3], targetY, Yeva, delta, denom, delta/denom, TMath::Power(delta/denom,2),fPenalty, reduYield); 
        //}
        if(kprint>3){
          exit(1);
        }
      }
    }
  }

  if(istep%1000==0){
    printf("istep %d myFCN %e pars %e %e %e %e %e -- %e %e %e\n", istep, myFCN(), localpar[0], localpar[1], localpar[2], localpar[3], localpar[4], localpar[5], localpar[6], localpar[7]);
  }
  istep++;

  return myFCN();
}

void MN2D::MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag)
{
  f = GetFCNChisquare(pars);
}

Int_t MN2D::DoMigrad(const Int_t nbeta, Int_t & counter, const Double_t tolDefault, const Int_t straDefault)
{
  //http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/node18.html   (List of Minuit commands)
  //ROOT's current policy is not to replace TMinuit with Minuit2
  //Minuit2 is not fully understood, like speed, error estimation.
  //TMinuit then need large tol, no need for minos. Anyway minos fails often.

  //let fMinuit exists outside if fit succeed, otherwise delete, return ierr
  //--> from Aug 4 2012 no check on convergence. The convergence is to varify smoothness behavior around the found minimum in a scale defined by tolerance, which can always be justified for convergence. The resulting minimum with/without convergence is no big difference. Actually smaller tolerance, even without convergence, has smaller minimum. With the use of decreasing regularization, the physical convergence is varified despite non-convergence in each fit.

  /*
http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_9.html
 MIGRAD

This is the best minimizer for nearly all functions. It is a variable-metric method with inexact line search, a stable metric updating scheme, and checks for positive-definiteness. It will run faster if you [SET STRategy]SET STRATEGY 0 and will be more reliable if you [SET STRategy]SET STRATEGY 2 (although the latter option may not help much). Its main weakness is that it depends heavily on knowledge of the first derivatives, and fails miserably if they are very inaccurate. If first derivatives are a problem, they can be calculated analytically inside FCN (see elsewhere in this writeup) or if this is not feasible, the user can try to improve the accuracy of Minuit's numerical approximation by adjusting values using the [SET EPSmachine]SET EPS and/or [SET STRategy]SET STRATEGY commands (see Floating Point Precision and [SET STRategy]SET STRATEGY). 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html
 SET STRategy<level>

Sets the strategy to be used in calculating first and second derivatives and in certain minimization methods. In general, low values of <level> mean fewer function calls and high values mean more reliable minimization. Currently allowed values are 0, 1 (default), and 2. 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_8.html
 MIGrad[maxcalls] [tolerance]

Causes minimization of the function by the method of Migrad, the most efficient and complete single method, recommended for general functions (see also MINImize). The minimization produces as a by-product the error matrix of the parameters, which is usually reliable unless warning messages are produced. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped even if it has not yet converged. The optional argument [tolerance] specifies required tolerance on the function value at the minimum. The default tolerance is 0.1, and the minimization will stop when the estimated vertical distance to the minimum (EDM) is less than 0.001*[tolerance]*UP (see [SET ERRordef]SET ERR). 

http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/node18.html
 IMProve 	[maxcalls]
If a previous minimization has converged, and the current values of the parameters therefore correspond to a local minimum of the function, this command requests a search for additional distinct local minima. The optional argument [maxcalls] specifies the (approximate) maximum number of function calls after which the calculation will be stopped. 

http://wwwasdoc.web.cern.ch/wwwasdoc/WWW/minuit/minmain/chapter2_9.html
 MINIMIZE

This is equivalent to [MIGrad]MIGRAD, except that if [MIGrad]MIGRAD fails, it reverts to [SIMplex]SIMPLEX and then calls [MIGrad]MIGRAD again. This is what the old [MIGrad]MIGRAD command used to do, but it was removed from the [MIGrad]MIGRAD command so that users would have a choice, and because it is seldom of any use to call [SIMplex]SIMPLEX when [MIGrad]MIGRAD has failed (there are of course exceptions). 
   */

  const Double_t tol = tolDefault;
  //const Int_t stra = straDefault;

  Double_t arglist[2];
  Int_t ierr = -999;

  Int_t stra = -999;
  for(stra = 2; stra>=0; stra--){
    printf("-----------> MN2D::DoMigrad counter %d tol %e loop in stra %d startnow ...\n", counter, tol, stra); 
    
    //always start from fMinuit=0x0, delete previous
    delete fMinuit;
    fMinuit = new TMinuit(Npar());
    
    fMinuit->SetFCN(MinFCNChisquare);
    fMinuit->SetPrintLevel(0);
    
    fMinuit->SetErrorDef(0.5);//ML 
    
    fStopwatch->Start(kTRUE);
    
    for(Int_t i=0;i<Npar();i++){
      const Double_t start=fPars[i];
      const Double_t ss=1e-2;
      const Double_t low= fLows[i];
      const Double_t hi= fHighs[i];
      
      fMinuit->DefineParameter(i, Form("p%d",i), start, ss, low,hi);
      if(fIfixed[i]){
        fMinuit->FixParameter(i);
      }
    }
    if(fMinuit->GetNumFreePars()!=Nfree()){
      printf("MN2D::DoMigrad fMinuit->GetNumFreePars()!=Nfree() %d %d\n", fMinuit->GetNumFreePars(), Nfree());exit(1);  
    }
    
    arglist[0]=stra;
    fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierr);
    
    arglist[0] = MaxIter();
    arglist[1] = tol;
    //fMinuit->mnexcm("MIGRAD", arglist, 2, ierr);
    
    //MINIMIZE is much better than MIGRAD
    fMinuit->mnexcm("MINIMIZE", arglist, 2, ierr);
    printf("-----------> MN2D::DoMigrad nbeta %03d counter %d tol %e finish ierr %d STRATEGY %d\n", nbeta, counter, tol, ierr, stra);
    
    if(ierr==0){
      break;
    }
  }

  //the results look quite similar to w/o IMPROVE  
  /*
  //to check real minimum
  if(fabs(fBeta-1)<EPSILON && ierr==0){
    arglist[0] = MaxIter();
    fMinuit->mnexcm("IMPROVE", arglist, 1, ierr);
    printf("-----------> MN2D::DoMigrad nbeta %03d counter %d tol %e finish ierr %d IMPROVE for beta=1\n", nbeta, counter, tol, ierr);
    //ierr=4 seems no new minimum, this is good
    if(ierr==4)
      ierr=0;
  }
  */

  printf("-----------> MN2D::DoMigrad nbeta %03d counter %d tol %e finish ierr %d finalstra %d\n", nbeta, counter, tol, ierr, stra);

  fStopwatch->Stop(); 
  fkErr = ierr;
  fStopwatch->Print("m"); //"m": milli sec; "u": micro sec

  //no check on convergence!
  /*
  if(ierr){
    printf("MN2D::DoMigrad fail %d %e\n", ierr, tol);
    delete fMinuit;
    //seems important!
    fMinuit = 0x0;
  }
  */

  counter++;

  //no modification to ierr
  //return ierr;
  return 0;
}

void MN2D::BetaFit(Int_t & nbeta)
{
  //do migrad with global fBeta
  //tol always start from IniTolerance
  //only modify fPars after loop

  if(fMinuit){
    printf("MN2D::BetaFit exit fMinuit already exist!!! It should only be initialized in DoMigrad!!\n"); exit(1);
  }
  
  Status(Form("-------> MN2D::BetaFit nbeta %d startnow", nbeta));
  
  Int_t mgcounter = 0;
  const Double_t tol=  dEdxGausN::IniTolerance();
  const Int_t stra = dEdxGausN::IniStrategy();

  //every time start from fPars which is not updated since entering BetaFit()
  DoMigrad(nbeta, mgcounter, tol, stra);
  
  Status(Form("-------> MN2D::BetaFit nbeta %d InLoopAfterDoMigrad mgcounter %d tol %e stra %d", nbeta, mgcounter-1, tol, stra));
  fStopwatch->Print("m");

  //now modify fPars
  //get par right after MIGRAD before HESSE modifies them
  Double_t dummy=-999;
  for(Int_t ipar=0;ipar<Npar();ipar++){
    fMinuit->GetParameter(ipar,fPars[ipar], dummy);
  }

  //update fStatChi2 and fPenalty with fPars
  GetFCNChisquare(fPars);

  Status(Form("-------> MN2D::BetaFit nbeta %d LastLoop mgcounter %d tol %e afterUpdatefPars", nbeta, mgcounter-1, tol/2.0));
  fStopwatch->Print("m");

  //save intermediate results and output to file!

  //important!!
  //only calculate the result for nbeta==0, i.e. no penalty!!
  //the error is affected by the penalty term, which reduce the effect of data statistics!!
  //i.e. errors of all nbeta should be the same as no penalty term!!
  //the error is determined by the contour of FCN: the larger the penalty term, the smaller the data affects the FCN, the smaller the error bars.
  //==> the error by the whole FCN describes the uncertainty of the line shape, the uncertainty of different NPolint, or different line shape should be consisdered as systematic error to be determined by comparing different sets of results
  //==> with the penalty defined as delta/StatErr, the total FCN has probability interpretation and the TMinuit err has correct meaning
  Int_t hesseflag = -999;
  Double_t hessearg[1]={MaxIter()};
  fMinuit->mnexcm("HESSE", hessearg ,1, hesseflag);
  printf("-------> MN2D::BetaFit nbeta %d hesse done. flag %d\n", nbeta, hesseflag);
  
  for(Int_t ipar=0;ipar<Npar();ipar++){
    Double_t tmperr = -999, dummy = -999;
    fMinuit->GetParameter(ipar, dummy, tmperr);
    fLowErrs[ipar]=tmperr;
    fHighErrs[ipar]=tmperr;
  }
  TMatrixDSym *covmat = BaseUtils::GetCov(fMinuit, fIfixed, fLowErrs);

  SetPhysics(covmat);

  Save(Form("outmn2d_post_iBeta%d", nbeta));

  nbeta++;

  delete covmat;
  delete fMinuit;
  //important: after delete, fMinuit can still have value. If later there is another delete, e.g. in Clear(), then it creshes!!
  fMinuit = 0x0;
}

void MN2D::Fit(const Int_t startnbeta)
{
  Status(Form("---> MN2D::Fit() startnow"));

  Int_t nbeta = startnbeta;

  fBeta = 0;
  //no penalty, should reproduce standard non-regularized results (with corresponding tol)
  //fStatChi2 and fPenalty are updated
  BetaFit(nbeta);

  //fNPolint!=0 : has penality
  if(fNPolint){
    if(fPenalty<EPSILON){
      printf("MN2D::Fit exit penalty 0 !! %e when fNPolint %d\n", fPenalty, fNPolint);exit(1);
    }
    
    //in case beta!=0, release all limits
    for(Int_t ipar=0; ipar<Npar(); ipar++){
      fLows[ipar] = 0;
      fHighs[ipar]= 0;
    }

    //with statistical error as the correct denominator in the penalty term, beta should just be trivially 1.

    const Double_t defactor = 2;

    // / 4. in kpenal makes it slower to compensate l_{reg} *0.5 due to complete chi2!! 17-march-2013
    Double_t kpenal = fPenalty;
    while(fStatChi2 > 0.5*kpenal){
      kpenal*=2;
      Status(Form("---> MN2D::Fit moderate penalty, increasing kpenal to %e", kpenal));
    }
    //very large penalty at first
    //start from chi2=beta*penalty ==> beta<1
    fBeta = TMath::Power(defactor, (Int_t)(TMath::Log(fStatChi2/kpenal)/TMath::Log(defactor)));
    Status((Char_t*)"---> MN2D::Fit very large penalty");

    while( fBeta<=1){
      if(nbeta>10)
        break;

      BetaFit(nbeta);

      Status(Form("---> MN2D::Fit InLoopAfterBetaFit nbeta %d", nbeta-1));
      fStopwatch->Print("m");
      
      fBeta *= defactor;
    }

    if(fkErr){
      fBeta = 1;
      Int_t normalFitN=100+startnbeta;
      
      BetaFit(normalFitN);
      Status(Form("---> MN2D::Fit finalBeta nbeta %d ierr = %d output as iBeta100", nbeta, fkErr));
    }
    else{
      fBeta /= defactor;
      if(fabs(fBeta-1)>EPSILON){
        printf("MN2D::Fit ierr strange fBeta %f\n", fBeta); exit(1);
      }

      Save("outmn2d_post_iBeta100");
      Status(Form("---> MN2D::Fit last Beta nbeta %d ierr = %d successful, copied as iBeta100", nbeta, fkErr));
    }
    fStopwatch->Print("m");
  }

  TFile *fout=new TFile(Form("outbeta%d.root", startnbeta),"recreate");
  fout->Save();
  fout->Close();
  delete fout;
}

